ON THE EXISTENCE OF ENERGY-PRESERVING SYMPLECTIC 
INTEGRATORS BASED UPON GAUSS COLLOCATION FORMULAE * 

LUIGI BRUGNANQt, FELICE lAVERNARO*, AND DONATO TRIGIANTE§ 

Abstract. We introduce a new family of symplectic integrators depending on a real parameter ol. For a = 0, the 
corresponding method in the family becomes the classical Gauss collocation formula of order 2s, where s denotes the 
number of the internal stages. For any given non-null a, the corresponding method remains symplectic and has order 
Is — 2: hence it may be interpreted as a 0(h^^~'^) (symplectic) perturbation of the Gauss method. Under suitable 
C " !) , assumptions, we show that the parameter a may be properly tuned, at each step of the integration procedure, so as 

^SJ ' to guarantee energy conservation in the numerical solution. The resulting method shares the same order 2s as the 

generating Gauss formula. 
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1. Introduction. Wc consider canonical Hamiltonian systems in the form 



(/ is the identity matrix of dimension m). Regarding its numerical integration, two main lines of 
investigation may be traced, having as objective the definition and the study of symplectic methods 
Cn ' and energy-conserving methods, respectively. In fact, symplecticity and the conservation of the 

^ , energy function are the most relevant features characterizing a Hamiltonian system. 

CO I From the very beginning of this research activity, high order symplectic formulae were already 

available within the class of Runge-Kutta methods, the Gauss collocation formulae being one notice- 
able example. One important implication of symplecticity of the discrete flow is the conservation of 

"O \ quadratic invariants. This circumstance makes the symplecticity property of a method particularly 

appealing in the numerical simulation of isolated mechanical systems in the form (|1.1|) . since it pro- 
vides a precise conservation of the total angular momentum during the time evolution of the state 
vector. As a further positive consequence, a symplectic method also conserves quadratic Hamiltonian 
functions (see the monographs [TH [THl [H] for a detailed analysis of symplectic methods) . 



, . , On the other hand, excluding the quadratic case, energy-conserving methods were initially not 

C^ ■ known within the class of classical methods and as a matter of facts, among the first attempts 

to address this issue, projection and symmetric projection techniques were coupled to classical 
non-conservative schemes in order to impose the numerical solution to lie in a proper manifold 
representing a first integral of the original system (see [131 Sect. VH.2], pQ [TT| and [T^ Sect. 
V.4.1]). 

A completely new approach is represented by discrete gradient methods which are based upon 
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the definition of a discrete counterpart of the gradient operator so that energy conservation of the 
numerical solution is guaranteed at each step and whatever the choice of the stepsize of integration 
(see 0110]). 

More recently, the conservation of energy has been approached by means of the definition of the 
discrete line integral, in a series of papers (such as [HI [16]), leading to the definition of Hamiltonian 
Boundary Value Methods (HBVMs) (see for example [3l|4]). These are a class of methods able to 
preserve, in the discrete solution, polynomial Hamiltonians of arbitrarily high degree (and hence, a 
practical conservation of any sufficiently difFerentiablc Hamiltonian J^. Such methods admit a Runge- 
Kutta formulation which reveals their close relationship with classical collocation formulae [5]. An 
infinity extension of HBVMs has also been proposed in [TUl [3] . 

Attempts to incorporate both symplecticity and energy conservation into the numerical method 
will clash with two non-existence results. The first |5] refers to non-integrable systems, that is 
systems that do not admit other independent first integrals different from the Hamiltonian function 
itself. According to the authors' words, it states that 

// [the method] is symplectic, and conserved H exactly, then it is the time advance 
map for the exact Hamiltonian system up to a reparametrization of time. 

The second negative result [7] refers to B-series symplectic methods applied to general (not neces- 
sarily non-integrable) Hamiltonian systems: 

The only symplectic method (as B-series) that conserves the Hamiltonian for arbi- 
trary H{y) is the exact flow of the differential equation. 

The aim of the present work is to devise methods of any high order that, in a sense that will 
be specified below and under suitable conditions, may share both features. More precisely, we will 
begin with introducing a family of one-step methods 

yi{a) = ^hiyo,a) (1.2) 

{h is the stepsize of integration), depending on a real parameter a, with the following specifics: 

1. for any fixed choice of a 7^ 0, the corresponding method is a symplectic Runge-Kutta method 
with s stages and of order 2s — 2; 

2. for a = one gets the Gauss collocation method (of order 2s); 

3. for any choice of j/o a-nd in a given range of the stepsize h, there exists a value of the 
parameter, say a*, depending on yo and h, such that H{yi) ~ H{yo) (energy conservation). 

As the parameter a ranges in a small interval centered at zero, the value of the numerical Hamiltonian 
function H{yi) will match H{y{tQ + h)), thus leading to energy conservation. This result, which will 
be formally proved in Section [H is formalized as follows: 

Under suitable assumptions, there exists a real sequence {ak} such that the numerical solution 
defined by yk+i = ^h{yk,ak), with yo defined in il.l]) . satisfies H{yk) = H{yo). 

To clarify this statement and how it relates to the above non-existence results, we emphasize that 
the energy conservation property only applies to the specific numerical orbit {yk} that the method 
generates, starting from the initial value j/o a-nd with stepsize h. For example, let us consider the very 
first step and assume the existence of a value a = ao, in order to enforce the energy conservation 



^We refer the reader to [2] for a complete documentation on HBVMs. 
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between the two state vectors j/q ^-^d ?/i, as indicated at item 3 above. If ao is maintained constant, 
the map y i-> $/j(j/, ao) is symplectic and, by definition, assures the energy conservation condition 
H{yi) ~ H{y()). However, it would fail to provide a conservation of the Hamiltonian function if 
we changed the initial condition yo or the stepsize h: in general, for any yo y^ yo, wc would obtain 
H{^h{yo, «o)) 7^ H(yo). Thus, the energy conservation property we are going to discuss weakens the 
standard energy conservation condition mentioned in the two non-existence results stated above and 
hence, by no means, the new methods are meant to produce a counterexample of these statements. 

The paper is organized as follows. In the next section we report the definition of the methods 
while in Section |3] we show their geometrical link with Gauss collocation formulae. In Section |4] we 
face the problem from a theoretical viewpoint and give some existence results that aim to explain 
the energy-preserving property of the new methods. In Section [5] we report a few tests that give a 
clear numerical evidence that a change in sign of the function g{a) = H{yi{a)) — -ff (j/o) does indeed 
occur along the integration procedure. 

2. Definition of the methods. Let ci < C2 < • • • < Cs and bi, . . . ,bshc the abscissae and the 
weights of the Gauss-Legendre quadrature formula in the interval [0, 1]. We consider the Legendre 
polynomials Pj{t), of degree j — I for j = 1, . . . , s, shifted and normalized in the interval [0, 1], that 
is 



Jo 



hj == 1, •••,«, 



(2.1) 



{Sij is the Kronecker symbol), and the matrix 



V 



( Pl(ci) P2(ci) 
PM) P2{C2) 



V Pl{Cs) P2{C 



Psici) \ 
Ps{c2) 



Ps{Cs) j 



(2.2) 



Our starting point is the following decomposition of the Butcher array A of the Gauss method of 
order 2s (see [131 Theorem 5.6]): 



A = VX,V-^, 



(2.3) 



where Xg is defined as 



X, = 



V 



L-i / 



(2.4) 



with 



2V(2j + l)(2j-l)' 



j = l,...,s- 1. 



(2.5) 
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We now consider the matrix Xs{a) obtained by perturbing (j2.4[) as follows: 



X4a) 



( \ -6 

V 



\ 



where a is a real parameter, and 



w. 



-(6-1+ a) 
6-1+ a 



/ 




X, + oiWs 



(2.6) 



-1 



V 



(2.7) 



1 y 

so that Xsipi) is a rank two perturbation of X^. 

The family of methods y\ = <f>/i(yo, o) we are interested in, is defined by the following tableau: 



(2.8) 



C\ 






A{a)=VXs{a)V~^ 


Cs 






bi bs 



Therefore 



and hence ^(0) ~ A. 



A{a) ^A + aVWsV 



(2.9) 



By exploiting Theorems 5.11 and 5.1 in [131 Chap. IV. 5], we readily deduce that the symmetric 
method (|2.8p has order 2s — 2 for any fixed a 7^ 0, and order 2s when a = 0. 



We set 



61 



bi 



n 



(2.10) 



bs J 



Theorem 2.1. For any value of a, the Runge-Kutta method defined in (|2.8p is symplectic. 

Proof On the basis of [HI Theorem 4.3, page 192], we will prove the following sufficient 
condition for symplecticity: 

nAia) +Aiafn^uju}'^. 



Since the degree of the integrand functions in ()2.1|) does not exceed 2s — 2, the orthogonality con- 
ditions may be equivalently posed in discrete form as 



^ bkPi{ck)Pj{ck) = Sij, i, j = 1, . . . , s. 



fe=i 
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or, in matrix notation, 



v'np = i. 



Considering that from (|2.11l) wc get V ^ = V'^fl, from (|2.9p we have that 

riAia) + Aiafn = nA + A'^n + anp{Ws + wDv'^^ = w ^"^ 

since the Gauss method is symplectic, and Wg is skew-symmetric so that Ws + W]" ~ 0. 



(2.11) 

(2.12) 
D 



In the event that a value a* 
conservation condition H{yi{a)) - 
method (|2.8p a symplectic scheme 



; a*{yo,h) for the parameter a may be found such that the 
H(jjo) be satisfied, we can extrapolate from the parametric 



y^<^h{y.a*), 



(2.13) 



that provides energy conservation if evaluated at yo. The existence of such an a* will be proved 
in Section H] One important implication the use of (|2.13p will guarantee is the conservation 
of all quadratic constant of motions associated with system (|l.ip . The fact that, in general, 
H{^ii{y,a*)) 7^ H{y), explains the extent to which the energy conservation property of the new 
formulae must be interpreted. Summarizing, the new formulae, when applied to the initial value 
system p.ip are able to define a numerical approximation of any high order, along which the Hamil- 
tonian function and all quadratic first integrals of the system are precisely conserved. 

2.1. Generalizations. The proof of Theorem I2.1l suggests how to extend the definition of the 
new formulae in order to get a family of methods depending on a set of parameters. Indeed, by 
looking at p.l2p . in order to preserve symplecticity, it is sufficient to substitute to the matrix aWs, 
any skew-symmetric matrix Ws of low rank, having non-null elements in the bottom-right corner so 
that, with respect to the Gauss method, the order is lowered as least as possible. 



Theorem 2.2. 



Consider the s x s matrix 



Vr 



where Vr is any skew- symmetric matrix of dimension r + 1 < s. The Runge-Kutta method defined 
by the Butcher tableau 



(2.14) 



Cl 


A{a) = V{Xs + Ws)V-^ 




bi bs 



is symplectic and has order p = 2(s — r). 

For example, a natural choice for the matrix Wg is: 

/ 

'■■ 



Ws 



a\ 



-ai 




(2.15) 



/ 
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leading to a multi-parametric method depending on the r parameters ai, . . . , a^- 

3. Quasi-collocation conditions. Condition (|2.9p reveals the relation between the Butcher 
arrays associated with the new parametric method and the Gauss collocation method. In order not 
to loose generality, just in this subsection we assume to solve the generic problem y ~ f{y)- 

We wander how the collocation conditions defining the Gauss methods are aff^ected by the 
presence of the parameter a. This is easily accomplished by expressing the coefficients of the 
perturbing matrix VWs'P~^ in terms of linear combinations of the integrals J' /j(r)dT, where Ij(t) 
is the jth Lagrange polynomial defined on the abscissae ci , . . . , Cs . Let F = (7^^ ) be the solution of 
the matrix linear system AT ^ VWsV^^, which means that (see p.Sp *) 

T = rX-^W,r-\ (3.1) 

The nonlinear system defining the block vector of the internal stages {Yi} is 

Y ^e(E>ya + h{A(E) I)F{Y) + ah{A F ® I)F{Y), 

where e is the vector defined in (j2.10p . hereafter / is the identity matrix of dimension 2m, and 

Y={Y,^ ... Y^)\ F{Y)^{f{Y,r ... f{Ys)^f- 

Therefore, the polynomial a{tQ + rh) of degree s that interpolates the stages Yi at the abscissae q, 
i = 1, . . . , s, is 

a{to + Th)=yo + hY, 'j(a;)dx/(r,)+a/i;^K]7fc, / lkix)dx] f{Y,). (3.2) 

j=i -^0 j=l \fe=i -^0 / 

Differentiating p.2p with respect to r gives 

s s / s \ 

a{to + rh)=Y, l,{r) f{a{to + c,h)) + « ^ E ^^M^) f^to + c,h))- (3.3) 

Finally, evaluating (|3.2I) at r = and p.3p at r = c^ yields 
CT(to) = 2/0, 



(j(<o + Cih) = f{a{to + Cih)) + a^jij f{a{to + Cjh)), z = l,...,s. 



(3.4) 



For a small, we can regard p.4p as quasi-collocation conditions, since for a = we recover the 
classical collocation conditions defining the Gauss method. 

3.1. Geometric interpretation. Let us assume the existence of a quadratic first integral 
M(y) independent from H[y): although this assumption is not strictly needed, it will somehow 
simplify the presentation of our argument. 

Roughly speaking, for a small, our parametric method may be interpreted as a symplectic 
perturbation of the Gauss method. Due to symplecticity of $^(., a), the parametric curve 

-^ = a^D^yi{a)^R^'^, (3.5) 
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a (L+xh) 

a ^ ' 

1 




H(y)=H(y, 



Fig. 3.1. A geometric interpretation of the parametric method 1)2. 8|l . Two quasi- collocation polynomials (Tai{io + 
rh) and ac,2{to +'rh) have as end-points the numerical solutions j/i(«i) and yi(«2) which are 0{h'^"~^) close to the 
Hamiltonian since, for a 7^ the method has order 2s — 2. This means that the length of the arc of curve enclosed by 
the points yi{ce\) and yi{a2) is 0{h^"~^). However, the parametric curve 7:06 [01,02] 1— > ?/i(q) passes through 
yi{0) which is at a distance 0(h^"~^^) from the manifold H{y) = H{yo), and there is a concrete possibility that this 
arc may intersect the manifold H{y) = H{ya) at a point yi(a*). 



where D is a given interval containing zero, will entirely lie in the manifold M{y) = M{yo) and its 
length will be 0{h'^''^^), sinee the method has order 2s— 2. However, the numerical solution produced 
by the Gauss method, namely 2/i(0) will be 0{h'^^~^^) close to the manifold H{y) = H{yo). Since the 
two manifolds contain the continuous solution, their intersection is nonempty and it is reasonable to 
expect that when a ranges in Z), yi{a) can slide from a region where iJ(j/i(a)) > Hijjo) to a region 
where H{yi{a)) < H{yQ), thus producing a sign change in the scalar function 



g{a)^H{y,{a))-Hiyo) 



(3.6) 



which, by continuity, will vanish at a point a*. 

Obviously, similar arguments can be repeated for the multi-parameter version (|2.14p of the 
method, where one has even more freedom in the choice of the parameters in the matrix Wg defined 
in (|2.15l) , in order to obtain the conservation of energyjj 

4. Theoretical existence results. After defining the error function (7(a) = H(jji{a)) — H{yQ), 
the nonlinear system, in the unknowns Yi, . . . , Ys and a, that is to be solved at each step for getting 
energy conservation, reads 



Y = e®yo + hiAia) ® I)F{Y), 
9{a) = 0, 



(4.1) 



■^We do not consider multi-parametric methods in the numerical results we present, since a single parameter 
suffices in getting the energy conservation property. 
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and its solvability is equivalent to the existence of the energy preserving method (I2.13|) we are looking 
for. After defining the vector function 

we see that system (|4.1[) is equivalent to G{h, yi,a) = 0. Of course G(0, yi,a) =0 for any value of 
a and, in particular G{0,yi,0) ~ 0. The Jacobian of G with respect to the two variables yi and a 
reads 

dG ., . f I 1^(2/0, a) \ 

-{h;yi,a) 



diyi,a) ' ' \^ V^H{yi) / ' 

where, as usual, / is the identity matrix of dimension 2m. From (|2.13|) we see that ^$^{yo,a) 
coincides with y'i{a) and, hence, with (T^(io + h). Due to the consistency of the method, it follows 
that, for a = 0, (7[^{tQ + h) -^ JVH{yo) as /i — s- 0. Therefore 

dG (I J"^H{yo) \ 

-(0,yi,0)= „_^, ^ :''' . (4.2) 



d{yi,a) \ V'^H{yo) 



Unfortunately, the Jacobian matrix (|4.2p is always singular. Consequently, the implicit function 
theorem (in its classical formulation) docs not help in retrieving existence results of the solution of 
(|4.ip when h is small. However, the rank of the matrix (j4.2p is 2m independently of the problem to 
be solved. This would suggest the use of the Lyapunov-Schmidt decomposition [22] that considers 
the restriction of the system to both the complement of the null space and the range of the Jacobian, 
to produce two systems to which the implicit function theorem applies. 

In our case this approach is simplified in that the implicit function theorem assures the existence 
of a solution Y{a) of the first system in (|4.ip for all values of the parameter a ranging in a closed 
interval containing the origin and \h\ < Hq, with Hq small enough. Then yi{a) = j/o + h{y^ (E) I)Y{a) 
is substituted into the second of (j4.ip to produce the so called bifurcation equation in the unknown 
a. When needed, we will explicitly write g{a,h) or _g(a, /i, j/o)i in place of g(a), to emphasize the 
dependence of the function g upon the stepsize h, that has to be treated as a parameter, and the 
state vector j/o- 

Let us fix a vector yo and look for solution curves of g{a, h) = in the (/i, a) plane. Obviously 
g{a, 0) = for any a, which means that the axis /i = is a solution curve of the bifurcation equation: 
of course, we are interested in the existence of a different solution curve a* = a*{h) passing through 
the origin. Since the gradient of g vanishes at (0,0), one has to compute the subsequent partial 
derivatives of g with respect to a and h. However one verifies that ij^ = ^^ = daOh evaluated at 
(0, 0) vanish as well, and this makes the computations even harder. For this reason, to address the 
question about the existence of a solution of (j4.ip . we make the following assumptions: 

(^i) the function g is analytical in a rectangle [—a, a] x [—h, h] centered at the origin; 
{A2) let d be the order of the error in the Hamiltonian function associated with the Gauss method 
applied to the given Hamiltonian system (|l.ip and the given state vector yo, that is: 

g(0, h) = H{y,{0)) - H{y„) = co/i'^ + 0{h'i+^), (4.3) 

with Co 7^ 0. Then, we assume that for any fixed a 7^ 0, 

g{a,h) =c{a)h'^-^+0{h'^-^), 



ENERGY-PRESERVING, SYMPLECTIC INTEGRATORS 9 

with c(a) ^ 0. 

Remark 4.1. A couple of quick comments are in order before continuing. Excluding the case 
where the Hamiltonian H{q,p) is quadratic (which would imply g[a,h) = for all alpha), the error 
in the numerical Hamiltonian function associated with the Gauss method is expected to behave as 
0{h^'^^^). Anyway, we cannot exclude a priori that special classes of problems or particular values 
for the state vector yo may occur, for which the order of convergence may be even higher. This is 
why we have introduced the integer d: therefore such integer will be at least 2s + 1. Moreover, we 
emphasize that the constant cq and the function c{a), will depend on yo. In conclusion, what we are 
assuming is that for the method (j2.8p , when a is a given nonzero constant, the order of the error 
H(yi{a)) — H{yQ) is lowered by two units with respect to the underlying Gauss method of order 2s, 
which is a quite natural requirement since such method has order 2s — 2. 

Theorem 4.2. Under the assumptions (Ai) and (A2), there exists a function a* = a*{h), 
defined in a neighborhood of the origin {—ho, ho), such that: 

(i) g{a*{h), h) ~ 0, for all h G {—ho, ho), 
(ii) a*{h) = const • /i^ + 0{h^). 

Proof From {Ai) and (^2) wc obtain that the expansion of g around (0,0) is: 

00 ^ r~.j 00 CX) ^ j-i^ -X- j 

^("-'^) = E-ya^(o^o)'^^ + E E 7i7ia^(0'0)'^^"- (4-4) 

We are now in the right position to apply the imphcit function theorem. Wc wiU look for a solution 
a* = a* {h) in the form a*{h) = r]{h)h'^, where rj^h) is a real-valued function of h. To this end, we 
consider the change of variable a = rjK^ , and insert it into (j4.4[) thus obtaining 



1 d'^q 

' ^7(0, Q)h'^^^ri + higher order terms. 



(4.5) 



[d-l)\dadh'i 
Therefore, for /i 7^ 0, g{a, /i) = is equivalent to g{r], h) = 0, where 



1 d'^g.^ ^. d' 



,d-l 



9 



1 d'^g 

+ - — — „, , , (0, Q)hri + higher order terms. 

d — 1 oaoh'^ ^ 



(4.6) 



By assumption (.42), both 077- (0, 0) and g^gj^d-i (0, 0) are different from zero and hence the implicit 
function theorem assures the existence of a function rj = ri{h) such that g{ri{h), h) = 0. The solution 
of g{oL, h) ~ for the variable a will then be given by 

a*(fe)=77(fe)fe^ = -— i— ^M^°'"^ h' + Oih'), (4.7) 

and this completes the proof. D 

By exploiting [T71 Theorem 6.1.2], we see that the function a*{h) is analytic if the power series 
is absolutely convergent for \h\ < ho and |a| < ao- In any event, the function a*{h) is tangent 
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to the h-axis at the origin which means that a very smaU correction of the Gauss method is needed 
when the stepsize is small enough. As a matter of fact, the needed correction is so small that the 
resulting method (|2.13p has indeed order 2s instead of 2s — 2, just as the Gauss method obtained 
by posing a = 0. This is a consequence of the following result. 



Theorem 4.3. Consider the parametric method p.Sp and suppose that the parameter a is 
actually a function of the stepsize h, in such a way that a{h) = 0{h^). Then, the resulting method 
has order 2s. 

Proof Let yi{a,h) be the solution computed by method (|2.8p at time to + h, starting at 
yo = y{to), and consider its expansion with respect to the variable a, in a neighborhood of zero: 

yi{a, h) = yi(0, h) + y'{C,a, h)a. 

We recall that yi{0,h) is the numerical solution provided after a single step of the Gauss method 
and hence it is 0(/i^*+^) accurate while, for a 7^ 0, yi{a,h) yields an approximation to the true 
solution of order 2s — 1. This implies that y'{Ca,h) is 0{h'^''^^). Consequently, 

yi{a, h) - y{to + h) = yi{0, h) - y{to + h) + y'{Ca,h)a = 0{h^'+^) + aO{h^'-^), 

from which we deduce that the error at the left hand side is 0(/i^''+^) if and only if a = 0{h^). D 

Figure WA\ reports the level curves of the function g(a,h) in a neighborhood of the origin, for 
the Kepler problem described in Subsection 15.11 (the vector yo has been chosen as in (|5.3p ) . The 
tick lines in the plot correspond to the points {a, h) in the plane where g vanish. This zero level set 
consists of the vertical axis h = and of the function a*{h), which splits the region surrounding the 
origin into two adjacent subregions where the function g has clearly opposite sign. Despite the local 
character of the above existence result, we see that the branches of the function a*{h) extend away 
from the origin. Similar bifurcation diagrams may be traced starting at different values of yo for all 
the test problems wc have considered: this suggests that, in the spirit of the long-time simulation of 
dynamical systems, a quite large stepsize may be used during the numerical integration performed 
by method (PTTO|) . 

We end this section by providing a straightforward generalization of Theorem 14.21 to the case 
where the parameter a is used to perturb a generic (not necessarily the last) element on the subdi- 
agonal of the matrix X^, and its symmetric. 

Theorem 4.4. Consider the method (|2.14p with Ws as in (|2.15p with a\= a and a2 = ■ ■ ■ = 
Ur = 0. We assume that assumption (Ai) and the following assumption (replacing (A2)) hold true: 

(A2) let d be the order of the error in the Hamiltonian function associated with the Causs method 
applied to the given Hamiltonian system (|1.1[) and the given state vector yo . That is, j^--^ 
holds true. Then, we assume that for any fixed a =^ 0, 



g{a, h) = c,(a)/i'^-2r ^ o(/,d-2r+i)^ 



with Cr{a) ^ 0. 



Then, there exists a function a* = a*{h) defined in a neighborhood of the origin (— /iq, ho) and such 
that: 

(i) g{a*{h),h) = 0, for all h G {—ho, ho), 
(li) a*{h) = const • h^"- + 0{h^''+^). 
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Fig. 4.1. Level curves in the plane {h, a) of the function g{a, h, j/o) associated with the method I I2.8II of order four, 
for the Kepler problem (see Subsection \5. l\ l . in a neighborhood of the origin: h G [—0.2,0, 2], a S [—0.5-10"'', 4-10~'']. 
Besides the a-axis, a zero level curve tangent to the h-axis at the origin is visible. Such curve separates two regions 
around the origin where the function g has opposite sign. We notice that just a small correction of the Gauss method 
suffices to recover the energy preservation even for relatively large stepsizes. 

The symplectic energy conserving method resulting from this choice of the parameter has order 2s. 
In the next section, we shall provide numerical evidence for the above presented results. 



5. Numerical tests. In this section we present a few numerical tests showing the effectiveness 
of our approach. Method (|2.13p and its generalization are implemented by solving, at each step, 
system (|4.ip . The efficient solution of such system will be the object of future studies; at present we 
adopt either one of the following techniques: 

1. at each step, an interval [011,02] is detected such that g{ai)g{a2) < 0; after that, a di- 
chotomic search is implemented to locate a* within an error close to the machine precision; 

2. the first (vector) equation in (|4.ip is solved with oq = (Gauss method) and oi = ch"^, 
where c and r are suitable constants empirically estimated^ after that, a sequence a^ is 
produced by solving the second (scalar) equation in (|4.ip via the secant method. 

In both cases, an outer iteration generating the sequence a^ converging to a* is coupled with an 
inner iteration that determines the solution yi{ak) starting from j/g. Such scheme is repeated at 
each step of integration. 



The methods that we will consider in our experiments are: method (j2.13p with s = 2 (fourth 
order); method (|2.13p with s = 3 (sixth order); the sixth-order method described in Theorem 14.41 
with s = 3, that is we insert a single perturbation parameter a in the first (rather than in the 
second) subdiagonal element of the matrix X^. In order to distinguish between these two methods 



•'For example see the last column in Table [5!T1 
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of order six, hereafter the latter wiU be referred to as "the order six method of the second type" . 

5.1. The Kepler problem. In this problem, two bodies subject to Newton's law of gravitation 
revolve about their center of mass, placed at the origin, in elliptic orbits in the (gi, g2)-plane. 
Assuming unitary masses and gravitational constant, the dynamics is described by the Hamiltonian 
function 

H{qi,q2,Pi,P2) = :^{pi+pi) , „ „ ■ (5.1) 

Besides the total energy H, a relevant first integral for the system is represented by the angular 
momentum 

L{qi,q2,Pi,P2) = qiP2-q2Pi- (5.2) 

Due to its symplecticity, the quadratic first integral (|5.2p will be automatically conserved by method 
(j2.8|) , for any choice of the parameter a. On the other hand, we show that, at each step of integration, 
the parameter a may be tuned in order to get energy conservation in the numerical solution. 

As initial condition we choose 



gi(0) = l-e, (Z2(0) = 0, pi(0)=0, ^2(0) = ^^^, (5.3) 

which confers an eccentricity equal to e on the orbit. Consequently, H{q,p) = —0.5 and L(q,p) ~ 
Vl — e^. We set e = 0.6 since, in this experiment, we are going to use constant stepsize (see [T^J 
Sec. 1.2.3]). More precisely, we solve problem (|5.ip in the interval [toj^] = [0,50] by the two-stages 
method (|2.13p with the following set of stepsizes: hi = 2^*, i = 1,...,7. Figure \EA] reports the 
errors in the Hamiltonian function H and in the angular momentum L of the numerical solutions 
generated by the method implemented with the intermediate stepsize h = 2~^. These plots, which 
remain almost the same whatever is the stepsize considered in the given range, testify that the 
integration procedure performed by method (|2.13|) is indeed feasible and both energy and angular 
momentum preservation may be recovered in the discrete approximation of (jl.ip . For comparison 
purposes, we also report the same quantities for the Gauss methods of order 4 (corresponding to the 
choice a = in (g^). 

The second and third columns of Table EH report the global error e{hi) = \yN{hi) — y{T)\, 
N — T/hi, at the end point of the integration interval and the corresponding numerical order. 
According to Theorem (|4.3p . we see that the maximum order is preserved by method (|2.13p . 



In Figure 15.21 the sequence a* , corresponding to the values of the parameter a that at each 
step restore the conservation of the energy, are plotted for the case h = 2~^. We consider d{h) = 
max„(a* ) — min„(a* ) as a measure of the total variability of the values of the sequence {a* }. Such 
quantity is reported in the fourth column of Table 15.11 for the values of the stepsize hi used in 
this test. According to the result of Theorem 14. 2[ the last column in the table confirms that the 
dependence of 5{h) on the stepsize h is of the form S = di? + h.o.t., with c ~ 0.16. 

5.2. Test problem 2. We consider the problem defined by the following polynomial Hamilto- 
nian function: 

H{qi,q2,Pi,P2) -= T^ipl+pl) + {ql + qlf- (5.4) 
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Fig. 5.1. Upper picture: errors in the Hamiltonian function of the Kepler problem evaluated along the numerical 
solution generated by the Gauss method of order four and its conservative variant (method 112.1311 with s = 2). Bottom 
plot: error in the numerical angular momentum of the solution computed by the two methods. In both cases the 
stepsize used is h = 2~^. 

This problem has been proposed in |19| as an example of a class of polynomial systems which, 
under suitable assumptions, admit an additional polynomial first integral F which is functionally 
independent from H. In this case, the additional (irreducible) first integral is 

L{qi,q2,Pi,P2) ^ qiP2- q2Pi- (5.5) 

The polynomial L being quadratic, we expect that our methods may preserve both H and Lo 



■*Of course L may again be interpreted as the angular momentum of a mechanical system having l|5.4|l as Hamil- 
tonian function. 
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Fig. 5.2. Sequence of the values of the parameter a* in the method ||2.13| I with s = 2 and h = 2 



h 


e(/i) 


order 


m 


S{h)/K' 


2-^ 


2.62 • 10" 




2.13-10-2 


8.5374-10-2 


2-2 


3.85-10-1 


2.763 


1.04-10-2 


1.6700 - 10-1 


2-3 


2.50-10-2 


3.945 


2.52-10-3 


1.6185-10-1 


2-4 


1.59- 10-3 


3.970 


6.23-10-'' 


1.5951-10-1 


2-5 


1.00 -10-'* 


3.991 


1.55-10-1 


1.5878-10-1 


2-6 


6.28-10-6 


3.997 


3.87-10-5 


1.5862-10-1 


2-7 


3.93-10-7 


3.999 


9.67-10-6 


1.5856-10-1 



Table 5.1 
Performance of the order four method l|2.13|l applied to the Kepler problem. The global error at T = 50 (second 
column), and the corresponding order obtained via the formula log2(e(/ii)/e(/ii-|_i)), indicate that the perturbations 
introduced in the Gauss collocation conditions (see l|3.4|l ) are small enough that the order 4 of the Gauss method with 
two stages is conserved by its energy preserving variant. The last two columns give a measure of the perturbations 
and of the rate they tend to zero as h —^ 0. The quantity S{h) is the amplitude of the minimum interval that encloses 
all the values a* for the given stepsize h and in the given integration interval. Hence the last column confirms what 
proved in Theorem\4.S\ namely that the perturbations are 0{h^). 



We have solved problem (|5.4p by means of two methods of order six (s — 3): method (|2.13p . 
and the order six method of the second type, described in Theorem 14.41 



Figure 15.31 reports the errors in the Hamiltonian function H and in the quadratic first integral 
L of the numerical solutions generated by the latter method implemented with the intermediate 
stepsize h = 2-3. For comparison purposes, we also report the same quantities for the Gauss 
methods of order six. 
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Fig. 5.3. Upper picture: errors in the Havailtonian function of test problem 2 evaluated along the numerical 
solution generated by the Gauss method of order six and its conservative variant of the second type. Bottom plot: 
error in the quadratic first integral I I5.5I I of the solution computed by the two methods. In both cases the stepsize used 
is h = 2-3. 



Tables 15.21 and 15.31 are the analogues of Table 15.11 for these two methods: we see that both 
methods achieve order six but, while in the former a*{h) = 0{h^), in the latter a*{h) = 0(h^) 
consistently with Theorems 14.21 14.31 and 14.41 



5.3. The Henon-Heiles problem. The Hcnon-Heiles equation originates from a problem in 
Celestial Mechanics describing the motion of a star under the action of a gravitational potential of a 
galaxy which is assumed timc-indcpcndcnt and with an axis of symmetry (the z-axis) (see |14j and 
references therein). The main question related to this model was to state the existence of a third 
first integral, beside the total energy and the angular momentum. By exploiting the symmetry of 
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h 


e{h) 


order 


S{h) 


S{h)/h' 


2-^ 


2.17- 10-^ 




1.59- 10-' 


6.37 -10-' 


2~2 


4.59-10-'^ 


5.562 


3.99-10-3 


6.39-10-2 


2-3 


7.77- 10-6 


5.884 


9.99-10-'' 


6.40- 10-2 


2-4 


1.24-10-^ 


5.970 


2.53-10--^ 


6.48- 10-2 


2-5 


1.94- 10-^ 


5.992 


6.33-10-5 


6.49- 10-2 


2-6 


3.05 -10-" 


5.994 


1.59-10-5 


6.51-10-2 



Table 5.2 
Performance of method 1 12.131 of order six applied to problerri 115.41 1. The reported quantities are the analogues 
of the ones presented in Table\5A\ 



h 


e{h) 


order 


S{h) 


S{h)/h'' 


2-^ 


4.91-10-2 




5.59-10-2 


0.895 


2-2 


1.46- 10-2 


1.753 


1.51-10-2 


3.87 


2-3 


1.84- 10-4 


6.304 


4.92 - 10-4 


2.01 


2-4 


3.23- 10-s 


5.836 


4.07 - 10-5 


2.66 


2-5 


4.73 -10-« 


6.091 


2.30-10-6 


2.41 


2-6 


7.03-10-1" 


6.074 


1.50-10-^ 


2.51 



Table 5.3 
Performance of the sixth-order method of the second kind applied to problem IJ5.4II . 



the system and the conservation of the angular inomentum, Hcnon and Heiles reduced from three 
(cyhndrical coordinates) to two (planar coordinates) the degrees of freedom, thus showing that the 
problem was equivalent to the study of the motion of a particle in a plane subject to an arbitrary 
potential U{qi,q2)'. 



Hiqi,q2,Pi,P2) = -^{Pl+Pl) + U{qi,q2). 



(5.6) 



In particular, for their experiments they chose 

U{qi,q2) = TiiQi + <l2) + <lil2 



^ll 



(5.7) 



which makes the Hamiltonian function a polynomial of degree three. When U{qi,q2) approaches 
the value h , the level curves of U tend to an equilateral triangle, whose vertices are saddle points of 

U (see Figure [Qjl. This vertices have coordinates Pi = (0, 1), P2 = (-^, -^) and P3 = (^, -i). 

Since U in (|5.6p has no symmetry in general, we cannot consider the angular momentum as 
an invariant anymore, so that the only known first integral is the total energy represented by ([5j 



ENERGY-PRESERVING, SYMPLECTIC INTEGRATORS 



17 




Fig. 5.4. Level curves of the potential U{qi,q2) of the Henon-Heiles problem (see Ij5.7|l ). The origin O is a stable 
equilibrium point, whose domain of stability contains the equilateral triangle having as vertices the saddle points Pi , 
P2, and Pi, provided that the total energy does not exceed the value g . Inside the triangle, a numerical trajectory 
(small dots) computed by the sixth-order method of the second type with stepsize h = 0.25 and in the time interval 
[0,500], is traced: its total energy is 0.15. 



itself, and the question is whether or not a second integral does exist. Henon and Heiles conducted 
a series of tests with the aim of giving a numerical evidence of the existence of such integral for 
moderate values of the energy H , and of the appearance of chaotic behavior when H becomes larger 
than a critical value: it is believed that for values of H in the interval (|, ^) this second first integral 
does not exist (see also [HI Section 13]). 



We consider the initial point Pq 
system a total energy H = 0.15 S (|, 



(gio, 920, Pio, P2o) = (0, 0, y ^, 0) which confers on the 
i ) . Therefore the orbit originating from Pg will never abandon 
the triangle for any value of the time t. We have integrated problem (j5.6|) in the time interval [0, 500] 
with stepsize h = 0.25 by using the Gauss method of order six and its conservative variant of the 
second type. Figure [53] shows the errors in the Hamiltonian function H in both cases. 



6. Conclusions. We have defined a new class of symmetric and symplectic one-step methods 
of any high order that, under somewhat weak assumptions, are capable to compute a numerical 
solution along which the Hamiltonian function is precisely conserved. This feature has been realized 
by first introducing a symplectic parametric perturbation of the Gauss method, and then by selecting 
the parameter, at each step of the integration procedure, in order to get energy conservation. A 
relevant implication of the symplectic nature of each formula is the conservation of all quadratic first 
integrals associated to the system. With the help of the implicit function theorem, we have shown 
that not only do these methods exist, but that the correction required on the Gauss method is so 
small that the order of convergence of this latter method is preserved by its conservative variant. A 
few test problems have been reported to confirm the theoretical results presented, and to show the 
effectiveness of the new formulae. 

This approach opens a number of interesting routes of investigation. First of all, if preferred, 
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Fig. 5.5. Errors in the Hamiltonian function | |5. 61 1- 115. 71 1 evaluated along the numerical solution generated by the 
Gauss method of order six and its conservative variant of the second type. Stepsize: h = 0.25; time interval: [0, 500]; 

initial condition (qio, (?20, PiO, P2o) = (0, 0, \/ ^, 0). 



the parameter could be selected in such a way to impose the conservation of other non quadratic 
first integrals different from the Hamiltonian function itself. More generally, the multi-parametric 
generalization introduced suggests the possibility of choosing the free parameters in order to impose 
the conservation of a number of functionally independent first integrals possessed by the continuous 
problem. Last but not least, the idea of considering symplectic corrections of the Gauss method 
could be in principle extended to other classes of symplectic methods known in the literature. The 
above described lines of investigation, as well as the efficient solution of the nonlinear systems arising 
from the conservation requirements, will be the subject of future researches. 
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